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ABSTRACT 

We study the formation of damped Lya and Lyman limit absorbers in a hierarchical 

clustering scenario using a gas dynamical simulation of an J7 = 1 cold dark matter 

universe. In the simulation, these high column density systems are associated with 

§^ forming galaxies. Damped Lya absorption, iVni ^ lO^"'^ cm~^, arises along lines of 

^ sight that pass near the centers of relatively massive, dense protogalaxies. Lyman 

O limit absorption, 10^^ cm~^ ^ Nm ^ lO^'''^ cm~^, develops on lines of sight that pass 

through the outer parts of such objects or near the centers of smaller protogalaxies. 

\^ The number of Lyman limit systems is less than observed, while the number of 

O damped Lja systems is quite close to the observed abundance. Damped absorbers are 

Q.^ typically ~ 10 kpc in radius, but the population has a large total cross section because 

O the systems are much more numerous than present day galaxies. Our results 

demonstrate that high column density systems like those observed arise naturally in 

a hierarchical theory of galaxy formation and that it is now possible to study these 

Qh absorbers directly from numerical simulations. 
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1. Introduction 

Damped Lya and Lyman limit absorption systems are among the most useful probes of highly 
nonlinear structure at moderate and high redshifts. They are much more common than quasars 
or luminous radio galaxies, and several lines of evidence suggest that they trace the population 
of "normal" galaxies in the young universe. The mass of atomic hydrogen in damped systems at 
2; ~ 3 is similar to the mass of stars in spiral disks today (Wolfe 1988), and observations show 
that at least two damped systems have linear extents ^ lOh^^ kpc (Briggs et al. 1989; Wolfe et 
al. 1993). Deep imaging and spectroscopic studies indicate that Lyman limit systems lie on lines 
of sight that pass near bright galaxies (Yanny 1990; Lanzetta & Bowen 1990; Bergeron & Boisse 
1991; Steidel, Dickinson & Persson 1994). Existing surveys have detected many Lyman limit 
absorbers (e.g., Sargent et al. 1989; Lanzetta 1991; Storrie-Lombardi et al. 1994; Stengler-Larrea 
et al. 1995) and damped Lya systems (e.g., Wolfe et al. 1986; Lanzetta et al. 1991), and the Keck 
telescope will push absorption surveys to even greater redshifts (Wolfe et al. 1995). 

Previous comparisons between observed properties of high column density absorbers and 
theories of structure formation have been based on simplified approximations. Mo & Miralda- 
Escude (1994) and Kauffman Sz Chariot (1994) use the Press- Schechter (1974) formalism to predict 
the abundance of dark matter halos and assume that all gas within these halos cools and becomes 
neutral. Ma & Bertschinger (1994) and Klypin et al. (1995) use a similar approach, but they 
normalize their Press-Schechter fits with dissipationless N-body simulations. These results can be 
compared to the amount of neutral gas inferred from the observations, but this analytic approach 
does not yield other properties of the absorption systems, and it can easily underestimate the 
constraining power of the observational data because it is based on the optimistic assumption of 
perfectly efficient gas cooling within the halo virial radius. 

In this paper we predict the absorption statistics of Lyman limit and damped Lya systems 
directly from a gas dynamical simulation, thereby eliminating many of the uncertainties required 
by previous methods. Our results for the Lya forest — lines with column densities below 
10^^ cm~^ — are presented in a companion paper (Hernquist et al. 1995; hereafter HKWM). 

2. The Simulation 

We explore a standard cold dark matter (CDM) universe with Q = 1, Hq = 50 km/sec/Mpc, 
fift = 0.05, and the perturbation spectrum normalized to crig = 0.7. We simulate a comoving 
periodic volume that is 22.22/(1 + z) Mpc on a side and include both radiative and Compton 
cooling for a gas of primordial abundance. Our initial conditions are identical to those of Katz, 
Hernquist & Weinberg (1992; hereafter KHW) and Hernquist, Katz & Weinberg (1995), but here 
the mass resolution is eight times higher: we use a total of 524,288 particles, half to represent 
the dark matter and half to represent the gas. The gas particle mass is 1.5 x IO^Mq and the 
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dark matter particle mass is 2.8 x 10^ Mq. Another important difference from KHW is that we 
include a photoionizing background, with uniform intensity J(^') = Jq {I'o/i') F{z), where uq is the 
Lyman-limit frequency, Jq = 10^^^ erg s^^ cm^^ sr^^ Hz~^, and 



We compute abundances assuming ionization equilibrium and an optically thin gas, and we use 
these abundances when computing rates of radiative cooling and photoionization heating. 

We originally performed this simulation to examine the impact of photoionization on the 
formation of low mass galaxies; the results of this study, and the details of the calculation, are 
described elsewhere (Weinberg, Hernquist Sz Katz 1995). Briefly, the simulation was performed 
with TreeSPH (Hernquist & Katz 1989), a Lagrangian code that combines smoothed particle 
hydrodynamics (SPH) with a hierarchical tree algorithm for computing gravitational forces. 
The cosmological version of TreeSPH and the methods used to compute radiative cooling in the 
presence of an ionizing background are described by Katz, Weinberg & Hernquist (1995, hereafter 



When computing gravitational forces, we use a critical opening angle, 6 = 0.7 {e.g. Barnes &; 
Hut 1986; Hernquist 1987). The equations of motion are integrated so that all the dark matter 
particles are moved on the system time step of 3.2 x 10^ years. TreeSPH, allows SPH particles 
to have time steps smaller than those of the dark matter particles by powers of two, up to a 
factor of 16 smaller than the system time step. The gravitational resolution is 20 comoving 
kpc (13 kpc equivalent Plummer softening), and the gas resolution varies from 5 comoving kpc 
in the highest density regions to 200 comoving kpc in the lowest density regions, reflecting the 
variable Lagrangian nature of TreeSPH. Compared to KHW, the higher mass resolution makes 
the simulation much more computationally expensive, so we have evolved it only to 2; = 2, taking 
300 CPU hours on a Cray C-90. 

When its power spectrum is normalized to crie = 0.7, the $7 = 1 CDM model reproduces 
the observed mass function of rich galaxy clusters reasonably well (White, Efstathiou & Prenk 
1993). However, this normalization is nearly a factor of two below that implied by the COBE 
microwave background fluctuations. This disagreement between the cluster-scale and COBE 
normalizations for "standard" CDM is, in our opinion, the most serious observational difficulty for 
this model. This discrepancy can be resolved in a number of ways, each being somewhat contrived 
but retaining a critical density universe dominated by cold dark matter. Among these variations 
are "tilting" or "breaking" the inflationary fluctuation spectrum, lowering the Hubble constant to 
Hq ~ 30 kms~^Mpc~^, including gravity wave corrections expected in typical inflationary models, 
or introducing a decaying particle that boosts the radiation content of the universe, delaying the 
epoch of matter domination. Each of these modifications would lower the amplitude of cluster-scale 
fluctuations relative to COBE fluctuations, and each would yield a somewhat different shape for 
the power spectrum on the scale of our simulation. Instead of selecting arbitrarily between these 
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possibilities, we opt to consider the standard CDM model with a sub-COBE normalization, so that 
we can compare our simulation to the wealth of earlier N-body and hydrodynamic studies of this 
scenario. We suspect that our results for the hydrogen absorption lines would translate without 
major modification to other = 1, CDM-dominated models that have a similar crie normalization, 
but the sensitivity of the hydrogen absorption lines to the shape of the power spectrum between 
cluster and sub-galactic scales is clearly an important issue for future investigations. 

To calculate the distribution in HI column densities of the absorption systems, we project 
the neutral hydrogen mass of each gas particle onto a 2-dimensional uniform grid using a cell 
spacing that is comparable to the highest resolution achieved anywhere in the simulation, about 
5.4/(1 + z) kpc. In this paper we consider absorption features with column densities in excess 
of 10^^ cm^^. The probability that more than one of them lies along a line of sight through the 
simulation is very small, since at 2; = 2, 3, and 4 the redshift interval across the simulation box, 
Az, is 0.019, 0.030, and 0.041, respectively, and, on average, there are observed to be only a few 
Lyman limit systems per unit redshift. To count systems with much smaller column densities one 
must construct artificial spectra, as in HKWM. 

Throughout the course of a simulation, we assume that the gas is optically thin to the ionizing 
radiation. This is adequate to compute the radiative cooling because gas that is dense enough 
to self-shield cools rapidly even when optically thin abundances are used. However, the neutral 
fractions are needed for the absorption line analysis, and self-shielding should be important in 
increasing the column densities of systems with Nhi > lO^^cm"^. In our analysis, we correct 
for self-shielding on a pixel-by-pixel basis, treating each absorber as a plane-parallel slab with a 
perpendicularly incident radiation field. While this is an idealized approximation, our conclusions 
are not sensitive to its details, which are described below. 

For each pixel we first calculate the neutral hydrogen column, neutral mass-weighted 
temperature, and neutral mass-weighted neutral fraction in the optically thin approximation. 
Using the temperature and neutral fraction, we solve for the average 3-dimensional hydrogen 
density using a bisection algorithm, then determine the total hydrogen column using the neutral 
mass-weighted neutral fraction and the neutral hydrogen column. To compute the self-shielding 
correction we assume that the total hydrogen is distributed uniformly in a plane-parallel slab, 
with half of the background radiation field incident perpendicular to each side of the slab. These 
assumptions do not themselves represent a realistic absorption geometry, but they are intended to 
be intermediate between a true slab and a sphere in an isotropic radiation background. 

We divide the slab into 100 equal computational bins. Starting at the outermost bin and 
sweeping through to the middle of the slab, we calculate new equilibrium abundances with the 
background field attenuated by HI, Hel, and Hell absorption in the other bins. However, since 
the abundances have now changed, the heating and cooling rates in the bin have also changed 
and hence there should be a new equilibrium temperature. We calculate the new equilibrium 
temperature including any adiabatic heating or cooling terms that are present. The new 
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temperature implies new equilibrium abundances, so we repeat the calculation, iterating within 
each bin until its abundances converge. We then repeat the computation for the full slab (since 
the new abundances imply new attenuation factors) and iterate until every bin in the slab has 
converged. 

At column densities of 10^*^ to 10^^ cm~^, the correction is small and varies between 1% 
and 10%. At higher column densities, the correction can be as large as a factor ^ 100. The 
self-shielding correction must be applied to over 200,000 pixels at each redshift, so we use the 
above procedure to create a 3-dimensional (column density, neutral fraction, temperature) lookup 
table and compute individual pixel corrections by interpolation. 

3. Results 

Figure 1 shows the projected neutral hydrogen map at z = 2. All lines of sight with column 
densities greater than 10^^ cm~^ appear white in Figure 1. These regions are well separated from 
one another. The interconnected, filamentary structures give rise to absorption at lower column 
densities. Even lines of sight through regions that appear black in Figure 1 can give rise to HI 
absorption at Lyo; forest column densities (see HKWM). 

We use the algorithm of Stadel et al. (1995: see also KHW, KWH) to identify gravitationally 
bound clumps of dense (p/p > 1000), cold (T < 30,000K) gas. These sites are likely to harbor 
rapid star formation and can therefore be identified with forming galaxies. They range from 
massive, high overdensity systems {p/p ^ 10^) containing hundreds of particles to lower overdensity 
clumps containing as few as 8 cold gas particles. The former are usually "mature" objects that 
began to form at 2; ~ 3 — 6, while the latter are young systems that have just begun to condense 
and cool. Figure 2 plots i^proj, the projected separation between an absorber and the center of the 
nearest galaxy in the simulation, against neutral column density, at z = 2. We sample a roughly 
equal number of lines of sight in each decade of -/Vhi. When I?proj exceeds 100 kpc, it is plotted at 
110 kpc. 

Figure 2 shows a clear association between high column density absorbers and forming 
galaxies. Damped Lya absorption, with TVhi > lO^'''^ cm~^, occurs along lines of sight that pass 
through the denser, more massive protogalaxies. The column density correlates (inversely) with 
projected separation, and the maximum -Dproj for damped absorption is about 20 kpc. Closer 
inspection of the neutral gas in these systems often shows evidence of flattening and rotational 
support, but our limited gravitational resolution of 20/(1 + 2;) kpc prevents us from drawing strong 
conclusions about the morphology of these objects. 

Lyman limit absorption (lO^'^cm"^ < iVni < 10^°'^cm~^) occurs on lines of sight that pass 
either through the outer parts (-Dproj ~ 20 — 100 kpc) of the more massive protogalaxies or near 
the centers of younger, lower density systems. Some of these absorbers have i^proj > 100 kpc, 
but these are mostly cases where the corresponding clump of cold gas has too few particles to 
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Fig. 1. — Projected neutral hydrogen in the simulation at z = 2. The simulation box is 22.22 
comoving Mpc across, which corresponds to 7.4 physical Mpc at this redshift. For the CDM 
cosmology used, this translates into 15.1 arc minutes. A self-shielding correction was applied to 
this HI map (see text for details). In the color scheme, white corresponds to TVhi ^ 10^^'^ cm^^, 
yellow to lO^^-^ cm-^ < Nm ^ lO^^-^ cm-^ red to lO^^-^ cm-^ < Nm ^ lO^^-^ cm-^ and black to 
Nm ^ IQi^-S cm-2. 
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be identified as a galaxy by our algorithm. If the simulation were repeated with better mass 
resolution, these lines of sight would probably have associated, low mass protogalaxies, and they 
might exhibit higher column density absorption because of denser gas and more efficient cooling. 

Figure 3 shows the HI column density distribution f{N, z), the number of absorbers per unit 
interval of column density per imit "absorption distance" AX = (1 + z)V2A2 (we use the Q = 1 
definition appropriate to our simulations instead of the = definition used in many observational 
papers). Histograms are computed from the projected neutral hydrogen maps at z = 2, 3, and 4 
by counting the fraction of pixels in each column density bin. The observational data for damped 
systems were kindly provided to us by K. Lanzetta (1995, private communication); earlier versions 
of these data were shown in slightly difi"erent form in Lanzetta (1991) and Wolfe et al. (1995). It is 
difficult to determine column densities for absorbers that are opaque in the Lyman continuum but 
not strong enough to produce damping wings. However, it is more straightforward to determine 
the cumulative number of systems with JVhi > 10^^'^, i.e. with an optical depth of unity at the 
Lyman limit. The diagonal boxes in Figure 3 represent the results of Storrie-Lombardi et al. 
(1994), which we have converted from cumulative numbers to constraints on f{N,z) by assuming 
that the form of f{N) is the iV~^'^^ power law found by Petitjean et al. (1993) in this regime of 
column density. The total redshift path length spanned by the 4096^ lines of sight through the 
simulation at the three different redshifts amounts to ~ 1.5 x 10*^, over 1000 times the total path 
length covered by all existing observations, though our lines of sight are not all independent. The 
simulated frequency distribution shows minimal evolution between z = 2 and z = 3 and weak 
evolution between z = 3 and z = 4:. 

4. Discussion 

For damped Lya absorbers, the agreement between the simulated and observed f{N,z) is 
fairly good. The simulation numbers are a bit low, but it is not clear that this difi'erence represents 
a significant failure of the CDM model. While our simulation volume is randomly chosen, it 
may not be large enough to be statistically representative — a larger box would include larger 
scale modes that would lead to the formation of "proto-Coma" clusters and "proto-Bootes" voids 
in which the collapse of galaxy scale objects would be, respectively, accelerated or retarded. 
Higher numerical resolution would probably increase the amount of cooled gas and hence f{N,z). 
However, star formation, neglected in this calculation, would reduce f{N,z) by converting neutral 
gas into stars, perhaps even reversing the sign of the evolutionary trend in Figure 3 (see, e.g., 
Kauffmann & Chariot 1994, Wolfe et al. 1995). We will examine these effects in future studies. 

Wolfe et al. (1995) estimate Qg, the density parameter of mass in damped Lya systems, 
by integrating f{N,z) between iVni = 10^°'^ cm~^ and A^hi = lO^^'^cm"^, the highest column 
density at which they have observed systems. They find Ug = 0.002 it 0.001 at z = 2, and 
^Ig = 0.004 lb 0.002 at 2 = 3, Since we know the physical state of all gas particles in the simulation, 
we can compute Qg directly by adding up the mass in cold, collapsed gas, i.e. gas with p/p > 1000 
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Fig. 2. — Projected distance i^proj between a line of sight with absorption column density iVni and 
the nearest galaxy in the simulation. Separations greater than 100 kpc are plotted as 110 kpc. 
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Fig. 3. — Logarithm of the frequency distribution f{N,z) vs. logarithm of the neutral hydrogen 
column density. Histograms show the simulation results at 2; = 2 (solid), z = 3 (dotted), and 2; = 4 
(dashed). Error crosses show observational results for damped sytems at z = 2 and z = 3, and 
diagonal boxes show observational constraints derived from the cumulative numbers of Lyman limit 
systems at 2; = 2, 3, and 4. 
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and T < 30, OOOK. We find Qg ^ 0.0065, 0.0036, and 0.0017 at z = 2, 3, and 4, respectively. 

At first glance it is surprising that our f{N,z) histograms lie below the observations but 
our values of are similar to or even larger than the Wolfe et al. (1995) values. However, since 
the f{N) distribution has a power- law index ~ —1.7, the contribution to Qg is dominated by the 
highest column density systems, and our directly measured Qg inchides gas with A^hi above the 
lO^^'^cm"^ upper bound used by Wolfe et al. (1995). The integration and extrapolation required 
to infer from the observations makes it difficult to disentangle observational/statistical issues 
from physical issues (the cosmological model, the effects of star formation). 

Wolfe (1988) and Schiano, Wolfe & Chang (1990) suggest that damped Lja absorption arises 
predominantly in large (~ 50 kpc) HI disks — the large size is needed to explain the observed 
incidence of absorption if the space density of damped systems is similar to that of present day 
disk galaxies. Our simulation reproduces the observed values of f{N, z) and Clg with objects that 
are only ~ 10-20 kpc in radius because the number of systems able to produce damped Lya 
absorption far exceeds the number of galaxies at ^; = 0. Some of these systems will likely 
evolve into dwarf galaxies, while others will merge to form the smaller number of larger galaxies 
present today. We shall be able to address this issue more quantitatively once we have evolved the 
simulation to z = 0. 

As shown in Figure 3, the number of Lyman limit absorbers falls short of the observed 
abundances by almost a factor of 10 at Nm 10^^ cm~^. Again, there is uncertainty caused 
by the small box size, but it seems likely that this discrepancy is statistically significant and 
could represent a failure of the CDM model we consider here. However, since many systems are 
marginally resolved, it may be that a modest increase in resolution would produce a large increase 
in the abundance of Lyman limit systems. We have compared this simulation to a run with the 
same initial conditions but 1/8 as many particles, and we find that the lower resolution has its 
most severe effect on the least massive radiatively cooled objects (Weinberg et al. 1995), which 
contribute much of the absorption near A^hi = lO^'^cm"^. Alternatively, we could be missing an 
additional population of absorbers produced by physical processes far below our resolution limit, 
such as fragmentation of cooling gas into pressure-confined clouds in galactic halos (e.g. Fall &; 
Rees 1985; Wang 1993; Mo 1994). 

A UV background field with Jq = 10^^^ erg s^^ cm^^ sr^^ Hz~^ is less intense than most 
estimates of the expected background from QSOs (e.g. Meiksin &; Madau 1993). We would 
not expect a stronger background to have much effect on the abundance of damped systems, 
since these would still be self-shielded and mainly neutral. However, a stronger background 
would substantially reduce the abundance of Lyman limit systems, worsening the confiict with 
observations. Conversely, increasing would increase the incidence of absorption by raising gas 
densities, neutral fractions, and cooling rates. 

Despite the many limitations of our numerical modeling, we have been able to show that the 
CDM model considered here reproduces important features of the observations: the column density 
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distribution (with the discrepancies already discussed) and the association between absorbers 
and galaxies. Equally important, we have demonstrated that it is now possible to predict the 
occurrence of hydrogen absorption directly from numerical simulations. This approach eliminates 
much of the uncertainty present in traditional semi-analytic calculations, so observations of 
damped Lya and Lyman limit systems can provide more robust constraints on cosmological 
theories. 

Although we have so far simulated only one cosmology, f2 = 1 CDM with aie = 0.7, we can 
draw at least one more general conclusion. A model having much less power on the scales of 
our simulation would produce even fewer absorption systems and would be in greater conflict 
with observations, especially if the formation of molecular gas or stars has depleted the atomic 
hydrogen. It seems likely, therefore, that any theory able to explain the amount of neutral gas in 
observed damped Lya systems must have more small scale power than crie = 0.7 CDM and/or 
a baryon density fifc > 0.05. This requirement may eliminate some theories that can otherwise 
provide a good account of galaxy clustering and fluctuations in the microwave background. 
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